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Chemical evolution of galaxies. I. A composition-dependent 
SPH model for chemical evolution and cooling 
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ABSTRACT 

We describe an SPH model for chemical enrichment and radiative cooling in cosmo- 
logical simulations of structure formation. This model includes: i) the delayed gas 
restitution from stars by means of a probabilistic approach designed to reduce the 
statistical noise and, hence, to allow for the study of the inner chemical structure of 
objects with moderately high numbers of particles; ii) the full dependence of metal 
production on the detailed chemical composition of stellar particles by using, for the 
first time in SPH codes, the Qy matrix formalism that relates each nucleosynthetic 
product to its sources; and iii) the full dependence of radiative cooling on the de- 
tailed chemical composition of gas particles, achieved through a fast algorithm using 
a new metallicity parameter C(T) that gives the weight of each element on the total 
cooling function. The resolution effects and the results obtained from this SPH chem- 
ical model have been tested by comparing its predictions in different problems with 
known theoretical solutions. We also present some preliminary results on the chem- 
ical properties of elliptical galaxies found in self-consistent cosmological simulations. 
Such simulations show that the above ^-cooling method is important to prevent an 
overestimation of the metallicity-dependent cooling rate, whereas the formalism 
is important to prevent a significant underestimation of the [cc/Fe] ratio in simulated 
galaxy-like objects. 

Key words: galaxies: abundances - galaxies: formation - galaxies: evolution - meth- 
ods: N-body simulations - methods: statistical. 



1 INTRODUCTION 

Chemical abundances are a rich reservoir of information 
which can help uncover how the formation of galaxies 
took place. Indeed, even after epochs where much in- 
formation is lost through phase mixing or violent relax- 
ation, stars still retain important clues about the evolu- 
tionary histories of the objects fro m which they come (e.g., 
iFreeman fc Bland-Hawthorn! [200^ ') . The earliest studies of 
the chemical evolution of galaxies were based on the wel l- 
known Closed Box model dLvnden-Bellll 19751 : iTinslevll 19801 ) . 
that played an important role in obtaining a first insight 
into this problem. This model considers galaxies as one- 
zone systems with a constant total mass and instantaneous 
recycling and mixing of matter. Because of such simplify- 
ing assumptions, the Closed Box model provides no idea 
about the internal structure of galaxies and leads to some 
predictions that are inconsistent with the distribution of G- 
dwarfs in the solar vicinity, as well as with other observations 
l|Pagel fc PatchettJI 19751 : (Valle et alj|20()3 h 



After the seminal work of lLacev fc Fall 

ll 19831 . Il985l~). other approa c hes (e.g. 
Il986l : iMatteucci fc Tornambel Il987l ; 



, .Diaz fc Tosi 
Clayton! Il987l 



ISommer-Larsen fc YoshiH Il989l : IMatteucci fc Francois! 
Il989l : iFerrini et al.lll994l : IChiappini et al.lll997l ) have con^ 
sidered multi-zone chemical evolution models with external 
gas infall. Numerical computations based on these kinds 
of models have been successful to describ e the radial 
distribution of abundances in our Galaxy (IFerrini et alj 
11994 iBoissier fc Prantzoj Il999h and in other spiral disks 
llMolla et al.lll996l : lBoissier fc Prantzodl200d : iMolla fc Dfaj 
l2005h . The success of most of these models is based on 
appropriate spatial variations of the ratio between the star 
formation rate and gas infall rate, which are achieved by 
using different code input parameters for different galaxy 
types. It is then important to compute more realistic models 
where such rates naturally appear, within a cosmological 
context, as a result of the physical processes involved in 
the formation and dynamical evolution of galaxies. Due to 
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the complexity of this problem, it must be addressed from 
approaches like, e.g., hydrodynamic N-body simulations. 

Among the different methods developed for the mod- 
elling of complex hydrodynamic phenom ena, the Smooth 
Particle Hydrodynamics (SPH) technique (Monaghan! ll992T l 
is one of the most widely used in astrophysics. Due 
to its Lagrangian nature, the evolving distribution of 
the gas, dark matter and stellar components can be 
easily followed in a self-consistent way. The first ap- 
proach to include chemical evolution in an SPH code 
posed bv ISteinmetz fc Mullen [1 1994 1 19951), 



by ISteinmetz fc Mullerl J 1994 119951). fol- 
lowed b vlRaiteri et al ] dl996h. ICarraro et al.l (1998). Bcrczik 
(|l999h . iKawata fc Gibson! (j2003h and iKobavashil |2004\ 
among others. All these works focussed on the chemi- 
cal enrichment of isolated objects, or groups of objects, 
formed fro m pre-prepared initia l conditions. Recent works 
(see, e.g., iMosconi et ail l200ll ; IScannapieco et aL 2005; 
iKobavashi et all l2007t ) have extended SPH simulations to 
study the detailed chemical enrichment of galaxies within 
a full cosmological context, where mergers and interactions 
have important effects. 

A new degree of sophistication in the SPH-modellin, 
of chemical enrichment was introduced bv iLia et al.l (|200 
hereafter LPC02), who proposed a stochastic algorithm to 
completely remove the assumption of instantaneous recy- 
cling of stellar ejecta. Such a delayed gas restitution from 
stars has non-negli gible dynamical effe c ts on N-body disks, 
as discussed by lJungwiert et al.l (|200ll . |2004| ) , and a direct 
impact on the resulting abundance gradients. The algo- 
rithm of LPC02 takes into account the non-instantaneous 
gas restitution from stars through a method where the to- 
tal number of baryonic particles remains constant and no 
dynamically hybrid particles are present: they are either 
fully collisional or fully collisionless. The presence of such 
hybrid particles would introduce spurious dynamical effects 
like, e.g., stars that follow for a while the evolution of the 
;as. The LPC02 approach, already used in some st udies 
Sommer-Larsen et al.1 [20051 : iRomeo et al.l [20051 . 120061 ). has 
the advantage that it leads to objects with reliable average 
values for their chemical properties. Nevertheless, due to its 
statistical nature, its main disadvantage is that large num- 
bers of particles per object are needed to avoid an excessively 
high scatter around the average values. 

Together with the inclusion of delayed gas restitution 
from stars, a necessary further improvement on the SPH 
modelling of chemical enrichment is the full dependence 
of both metal production and radiative cooling on the de- 
tailed chemical composition of star and gas particles, re- 
spectively. Indeed, the metal production of stars is usually 
included in SPH codes on the base of up-to-date libraries 
with stellar models and their corresponding ejecta for stars 
with different masses and metallicities. S uch libraries (e.g., 
Portinari et al.|[l998l ; [Gavilan et al.ll2005l ) generally assume 
solar relative abundances of the various species for any total 
metallicity Z. Nevertheless, it is well known that within a 
given Z different chemical compositions are possible. Chem- 
ical evolution models show in fact that abundance ratios 
are not constant in the c ourse of galactic evolution (see, 
e.g., IPortinari et alj Il998l ). Deviations from solar propor- 
tions, particularly the ratios of a-peak to Fe-peak elements, 
are also observed in different stellar systems. Supra-solar 
values for these ratios are observed in metal-poor stars of 



our Galaxy as well as in g lobular clusters l|Kuntschner et all 
120021 : lLarsen et al.l |2002| ) , where [a/Fej seem s to increase 
with decreasing metallicity l|Puzia et al.ll2005l) . reaching val- 
ues that vary from some a elements to others: Mg seems 
to conver ge towards [Mg/ Fe] ~ 0.4 for very metal-poor 
stars (see lFuhrmarm1ll999|) wh ereas oxygen trends towards 
[O/Fe] ~ 0.6 (see lReetzlll999l ). Such an a-enhancement is 
also observed in spiral galaxy bulges (jSansom et al.|[l994l ) 
as well as in most of the c entra l regions of ell i ptical 
galaxies llBurstein et ail 1 1984 1 19881: IWorthev et al.l Il992l: 



Guzman et al.TTl99 



Carollo et al.l 1 19931 : iDavidge fc Clarkl 
1994 IWorthevlll994 ). where the [a/Fe] ratio increases with 
velocity dispersion and hence with mass. Since the a- 
cnhancement is commonly interpreted as the result of large 
and short star formation events at early times, it could pro- 
vide us with a very strong constraint for models of galaxy 
formation. To accurately follow the evolution of this or any 
other abundance ratio, the assumption of solar proportions 
must be relaxed in cosmological simulations when comput- 
ing the production of the different chemical elements. 

On the other hand, the chemical content of the gas in 
the interstellar medium (ISM) is enriched by metal ejecta 
of different stars and then mixed through complex processes 
involving both local diffusion and motions at larger scales. 
Consequently, the gas of cosmological simulations (as well 
as in real objects) has a non-trivial mixture of metals with 
trends that are not necessarily the same as in the solar vicin- 
ity. Since the radiative cooling rate of gas closely depends 
on its metal composition, a realistic modelling of galaxy 
formation must consider a metal-dependent cooling func- 
tion. As pointed out by some authors (e.g., I Romeo et al.l 
I2006D . ideally one would need cooling functions for different 
relative proportions of different elements. Nevertheless, the 
lack of a fast algorithm to implement such a composition- 
dependent cooling rate has forced the u se in c osmological 
simulations of the ISutherland fc Dopital (|l993l ) functions, 
given for chemical compositions mimicking those in the so- 
lar vicinity (scaled, according to the total metallicity, from 
primordial to solar abundances). 

In this paper we present a new SPH-model for chemical 
evolution and cooling. Our model includes: i) the delayed gas 
restitution from stars by means of a probabilistic approach 
based on the LPC02 scheme. Such a scheme has been modi- 
fied to reduce the statistical noise and, hence, to allow for the 
study of the inner chemical structure of objects with mod- 
erately high numbers of particles; ii) the full dependence of 
metal production on the detailed chemical composition of 
stellar particles. To this end, the chemical production of a 
stell ar particle is computed b y using the matrix formal- 
ism (|Talbot fc Arnettlll973l ). that relates each nucleosyn- 
thetic product to its sources; and iii) the full dependence 
of radiative cooling on the detailed chemical composition of 
gas particles. This latter issue is achieved through a fast al- 
gorithm where the cooling rate is not computed by using the 
total metallicity as the scaling parameter. We use instead a 
parameter C(T) defined as a linear combination of the abun- 
dance of different chemical species. The coefficients of such 
a linear combination depend on temperature and give the 
weight of each element on the total cooling function. 

The work presented in this paper will focus on the de- 
scription and testing of our chemical model. In order to 
analyse whether this model is able to produce galactic ob- 
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jects with chemical properties consistent with observations, 
we also report some first results on the cosmological forma- 
tion of elliptical galaxies. Such simulations have b een carried 
out by using the DEVA code |Serna et al.ll2003h . Neverthe- 
less, the model presented in this paper is not limited to any 
particular code. In forthcoming papers we will analyse for 
larger galaxy samples some important chemical properties 
like, e.g., the metallicity distribution functions (MDF) and 
abundance gradients within individual objects, as well as 
the mass (or luminosity) relation with age, metallicity and 
[a/Fe]. 

This paper is organised as follows. In Section 2 we de- 
scribe our model for chemical evolution, whereas Section 3 
presents our algorithm to compute composition-dependent 
cooling functions. Some synthetic tests as well as first results 
on cosmological simulations are shown in Section 4. Finally, 
in Section 5 we summarise our main results. 



2 THE CHEMICAL EVOLUTION MODEL 

As in any particle-based scheme, both the dark matter and 
baryonic content of the simulation box are sampled in our 
model by using a discrete number TV of particles. Each bary- 
onic particle can be either in the form of gas or in the form 
of a stellar particle representing a single stellar population 
(SSP) of mass m and age t. We do not consider hybrid par- 
ticles simultaneously hosting both a stellar and gas content. 
Star formation (SF) and gas restitution from stars are then 
modelled by turning gas into stars, or stars into gas, respec- 
tively. 

Within the above scheme the chemical evolution model 
must provide us with procedures to compute: i) when and 
how gas particles are turned into stars (star formation); ii) 
when and how stellar particles are turned into gas (gas resti- 
tution); iii) the metal production of each stellar particle over 
a timestep At (metal production); and iv) how the chem- 
ical production of stars is released and mixed through the 
gas component (metal ejection and diffusion). We will now 
address each of these issues separately. 



2.1 Star Formation 



Star formation is 
codes for galaxy 
been proposed in 
how g as particles 



commonly included in most SPH 
formation. Different criteria have 
the literature to decide when and 
turned into stars (e.g ., iKatzl 



19921; ISteinmetz fc Mulled Il994l; lYepes et al 



Berczik fc Petrovl l200ll; ISpringel fc Hernquistl 



1997 



2003 



Merlin fc Chiosil l2007f ). In our current implementation 



of SF we consider that gas particles are eligible to form 
stars if they are located in a region with a convergent flow 
and a gas density higher than a given threshold, p t h- In 
that case they are turned into stellar particles according 
to a Kennicutt-Schmidt law-like transformation rule (see 
lKennicutdll99i Isiikll200ll ). 



dpg 

dt 



dp* 
~dt 



c *Pg 

to 



(1) 



where p g and p t are gas and stellar density, respectively, c, is 
a dimensionless star-formation efficiency parameter, and t g 



is a characteristic time-scale chosen to be equal to the max- 
imum of the local gas-dynamical time td yn = (47rGp g ) _1//2 , 
and the local cooling time, t coo i — Ui/iii, where Ui is the 
internal energy. Eq. Q implies that the probability p that 
a gas particle forms stars in a time At is 



i — c*At/t„ 

p = 1 — e 19 



(2) 



As usual, we compute p at each time step for all eligible 
gas particles and draw random numbers to decide which of 
them form stars in the time interval [t,t + At]. Then, each 
of these randomly selected gas particles is transformed into 
one stellar particle. 



2.2 Gas restitution 

Once stars are present they return to the ISM part of their 
mass in form of chemically processed gas. As already men- 
tioned, each stellar particle is treated as an SSP of total 
mass m and age t. Within each SSP, stellar masses are dis- 
tributed according to a given Initial Mass Function (IMF) , 
<fr(M). Throughout this paper we use the lChabrie"rl (|2003ah 
IMF with a mass range of [Mi, M u ] = [0.1, 100] M w . Such 
an IMF is similar t o other possible choices (e.g., ISalpeten 
1 19551 ; iKroupal fl998l) . but provides a better fit t o counts of 
low-mass stars (|Chabrieill2003b1 ; iBell et al.ll2003h . 

Individual stars of mass M are characterised by a mean- 
lifetime r(M). Therefore, within an SSP of age t, any star 
with t(M) < t has already died so that part of its mass 
remains as a stellar remnant, Mr, while the rest should be 
in the form of gas ejected back to the ISM. The gas mass 
fraction, E(t), of an SSP of age t is then given by: 



E(t) = 



M(t) 



M - M r {M) 
M 



$(M)dM 



(3) 



where M u is the upper mass limit for stars in the IMF, and 
M(t) is the mass of stars with lifetime t so that any in- 
dividual star with M > M(t) has already ejected a mass 
M — M r (M) of gas. We compute ((3J by using the stellar 
mean-lifetimes r(M) from the Geneva evolutionary tracks 
(Schallcr et al. 1992) and the remnant masses M r (M) of 
iGavilan et aTTi|2005h for M < 8M Q , or IWooslev fc Weaver! 
<|l995h for M > 8M Q . 

In a stochastic approach the above mass fraction of gas 
can be used to compute, at each timestep At, the probability 
that an entire stellar particle turns back into gas through a 
Monte Carlo method similar to that used for SF. Such a 
probability is given by (LPC02): 



E(t + At)-E(t) _ f* +At e(t')dt' 



Pa 



1 - E(t) 



1 - E(t) 



where 



e(t) = 



dE(t) 
dt 



(4) 



(5) 



is the instantaneous rate of total ejecta, while [1 — i?(t)] -1 
is a correction factor that accounts for the stellar particles 
that have already turned to gas, so that the probability p g 
is not computed for them anymore. 
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2.3 Metal production 

We assume that initially all gas particles have a primordial 
composition. As part of stellar evolution, newly produced 
elements are released to the surrounding ISM, either as stel- 
lar winds or byproducts of supernovae (SNae) explosions. 
We consider the evolution of the following elements: H, 4 He, 
12 C, 13 C, 14 N, ie O, 20 Ne, 24 Mg, 28 Si, 32 S, 40 Ca and 56 Fe. 

For an SSP of age t, the instantaneous ejection rate of 
a given element i can be expressed as: 



i(t) = Pi (t) + e(t)Xi 



(6) 



where pt(t) is the ejection rate of the newly synthesised ele- 
ment i (i.e., the yield of i), while e(t)Xi is the ejection rate 
of the element i already present when the SSP was formed 
from a gas particle with abundance Xi. All these rates are 
expressed as quantities per total mass unit. 

As already quoted in § [TJ the yield of any element i 
depends on the detailed abundances of all other species j. 
In order to take into accou nt such a dependence we use the 
Qij formalism proposed bv lTalbot fc Arnettl (1 19731 ). Such a 
formalism links any ejected species to all its different nu- 
cleosynthetic sources, allowing the model to scale the ejecta 
with respect to the detailed initial composition of a star. 

For an individual star of mass M, the mass fraction 
initially in the form of chemical species j, transformed and 
ejected as chemical species i, is written in the Qij formalism 
as: 



Q,(M) = *g= 



Mi. 



XjM 



(7) 



where Xj is the initial abundance of j, whereas Mij,e Xp is the 
mass of i that was synthesised from j and finally expelled. 

For a whole SSP, the contribution to ei(t) due to j is 
then given by 



©,•(*) = -Q«(M(t))*(M(t))^ 



(8) 



and, therefore, the total ejection rate of the chemical species 
i can be written as: 



0) 



and the yield of i can be obtained from Eqs. [6] and [9] as 

Pi(i) = e<(t) - e(t)Xi . (10) 

To compute qij(t) we consider the element production 
from: 

(i) Enrichment from low and intermediate mass (LIM) 
stars. We assume that all stars with masses in the range 
0.8 - 8 M Q end their life, after an AGB or TP-AGB phase, 
by the loss of their envelope that is ejected to the ISM as 
enriched gas in the for m of planetary nebu lae. We use the 
set of stellar yields bv IGavilan et a l. (2005), which include 
the effects of the third dredge-up for the TP-AGB stage 
and of the hot-bottom burning processes. In this set the pri- 
mary nitroge n contribution of the LIM stars is higher for low 
metallicities (|Gavilan et alj [20061 ) and, simultaneously, the 
total nitrogen yield is lower than in o ther ste l lar yi eld sets 
ijvan den Hoek fc Groeneweeenl 1 19971 : iMarigol l200lh . Mod- 
els using these stellar yields produce results in excellent 
agreement with observation of nitrogen from our Galaxy 



(halo and disk stars), from extragalac tic HII regions, and 
from Damped-Ly man- Alpha galaxies (jGavilan et al.1 120061 : 
iMolla et al.1l2006i ). 

(ii) Type II supernovae (SNII). We assume that stars 
more massive than 8 Mq produce S NII. For the chemi- 
cal pr oduction we adopt the yields of IWooslev fc Weaver! 
l|l995l 'Fl These yields include the elements produced in pre- 
supernova evolution and SNII explosions for metallicities be- 
tween Z — and Zq a nd masses between 11 and 40 M(rjf|. 
Unlike other sets (e.g., iPortinari et al. I ll998l ). these yields 
do not consider t h e con t ributi on of stellar winds. As shown 
by IGavilan et al. I (|2005l . I 2006), when stars suffer mass loss 
from winds they eject large quantities of helium and carbon 
from their envelopes. This hinders the creation of oxygen 
and other heavier elements, and results in a low yield of 
oxygen and a high yield of carbon. 

(iii) Type la supernovae ( SNIa). We use the c hemical 
yields from the W7 model bv llwamoto et all (|l999l ). tabu- 
lated for two different metallicities (sol ar and sub-solar) . The 
SNIa rates are computed according to iRuiz-Lapuente et al.l 
( 2000) who provided us with a numerical table (private com- 
munication) with the time evolution of the supernova rates 
for a single stellar population. These rates were derived for 
several combinations of possible candidates of binary system 
or SNIa (double degenerate, single degenerate, etc.), consid- 
ering all parameters (secondary lifetimes, orbital velocities, 
distances between both stellar components, etc.) that deter- 
mine the conversion of a binary system into an SNIa explo- 



All the above sets of stellar yields are used as input data 
to construct the gy(t) matrix a ccording to thelFerrini et al.l 
||1992| ) algorithm as updated in IGavilan et ail (|2005|). Such 
an updated algorithm is based on iGalli et al.l (|l995f ) and 
IPortinari et all \ 19981) . 



2.4 Metal ejection and diffusion 

2.4-1 Metal ejection 

How the chemical elements are distributed and mixed in 
the ISM is a complex problem that takes place at scales 
much smaller than the resolution reached in cosmological 
simulations. Most SPH implementations of chemical enrich- 
ment do not consider gas restitution from stars. Therefore, 
at each timestep At, the metals produced by a stellar parti- 
cle are distributed through the neighbouring gas by means 
of, e.g., SPH-spreading. This simple procedure fails however 
when the gas restitution from stars is considered. Indeed, 
the most straightforward way of restoring gas is to perform 
the opposite process of star formation, and turn whole stel- 
lar particles into gas. However, when a stellar particle of age 
t turns back into a gas particle, its metal production from 
stars with lifetimes longer than t has not been considered 
yet. If, at that timestep, such a remaining metal production 
is computed and distributed over the neighbouring gas, par- 
ticles with an extremely high metallicity could result. Oth- 



1 The Fe ejecta has been divided by 2, following T immes et al.l 
lll995l) 

2 We extrapolate linearly for ejecta up to the upper mass limit 
considered (100 M©). 
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erwise (i.e., if the remaining metal production is ignored), 
simulations including gas restitution would systematically 
underestimate the chemical enrichment. A straightforward 
procedure to circumvent this problem consists in introducing 
new gas particles and/or modifying the individual masses of 
particles to progressively account for the gas restitution from 
stars. 

Within a scheme where both the individual mass and 
the total number of particles remain constant, as we con- 
sider here, the metal release from stars can be also addressed 
through a statistical approach like that proposed by LPC02. 
That is to say, we must evaluate the overall metal produc- 
tion of a large enough number N of stellar particles with 
similar ages and initial metallicities, so that they themselves 
constitute an SSP. The way in which individual particles re- 
lease their metals must then ensure the right overall chemical 
evolution for the whole set of N particles. For example, a 
possible procedure suggested by LPC02 is based on the fact 
that, within a timestep At, a stellar particle of age t theoret- 
ically ejects an amount ' e(t')dt' of gas with chemical 
abundances given by 



AX'i(t) = 



This procedure then considers that, as long as a parti- 
cle remains as a stellar one, there is no gas restitution and 
metals are not released to the surrounding gas. On the con- 
trary, if such a particle turns back into gas at some timestep 
At, the chemical content of the resulting gas-again particle 
is updated so that it corresponds to the composition of the 
gas released at such timestep, given by Eq. pip . Over a suf- 
ficiently large number N of particles, constituting a whole 
SSP, this approach gives a fair representation of the over- 
all metal production. The main limit of this method is that 
newly produced metals for a whole set of stellar particles are 
sampled by the few particles of different ages that turn back 
into gas within a given timestep, i.e., within a given cosmic 
age interval, [t, t+ At]. Since the chemical composition of the 
gas released from an SSP is a function that strongly varies 
with time (see Fig. [TJ , a large number of particles is then 
needed to obtain a meaningful sampling, but also to avoid 
strong statistical noise which could lead to excessively high 
scatter in the resulting distribution of abundances. 

In order to reduce the above statistical noise we have 
implemented a different procedure, in such a way that the 
metals produced by any stellar particle at a given timestep 
[t, t + At] are distributed through the neighbouring gas at 
this timestep (this procedure is similar to that used in pre- 
vious SPH implementations of chemical enrichment with- 
out restitution). Actually, since we incorporate the diffusion 
of metals through the gas (see below), just the nearest gas 
neighbour receives the yields of a stellar particle during each 
timestep. Note that, since stars do not follow the motion of 
gas particles, the identity of such nearest gas neighbour can 
change from one step to another. The composition X[ of 
the nearest gas neighboui[f| of a stellar particle of age t is 
updated according to: 

3 When a stellar particle is turned into gas, the resulting gas- 
again particle is considered as the nearest gas neighbour in that 
timestep. 



J t t+at Mt')dt' 

1-E(t) 



(12) 



Here the yield pi (t) is given by f| 10|) and the factor [1 — 
Z?(t)] _1 converts yields per total mass unit to yields per 
stellar mass unit and then corrects for the missing metal 
production of stellar particles that were already transformed 
into gas. To see that is indeed the case, consider No particles 
of mass m that become stellar particles at to = 0. Such 
No particles can be regarded as an SSP. For each chemical 
species i, the mass production by this SSP during the time 
interval [t, t + At] is: 



AM, (t) =mN(t)AX'i(t) , 



(13) 



where mN(t) — mNo[l — E{t)\ is the actual stellar mass 
available at time t, after a fraction E(t) of the stellar mass 
initially available at t = has been transformed again into 
gas. Taking into account (|12|1 . we can write: 



AMj(t) = mN I pi(t')dt' 



(14) 



and this is just the expression for the theoretical mass pro- 
duction during the time interval [t, t+ At] of the i-th element 
by the whole stellar mass in the SSP, mNo. 

Note also that, since the yield definition implies 
EiPi(*) = 0, Eq. Gl gives: 



^AM,(i)=0 



(15) 



Therefore, when the composition of a gas particle is up- 
dated using (|12|) . there is no mass exchange between such 
a gas particle and the neighbouring stellar particle that re- 
leases its metals. The mass of each particle remains then un- 
altered and the only mechanism transforming stellar mass 
into gaseous mass is still that described in § !2.2l for gas resti- 
tution. 



2.4-2 Metal diffusion in the gas 

Once the stellar production of metals has been released to 
the ISM, their redistribution through the gas component is 
governed by the turbulent motions of the gas. The turbulent 
mixing that takes place at subresolution scales can b e prop- 
erly modelled by a diffusi on law at resolved scales (jTavlorl 
Il92ll : iKlessen fc Lmll2003l ) . 

^ = DV 2 X t . (16) 
at 

The SPH formulation of the d iffusion equation for a 
compressible fluid has been given by iMonaghanl (|2005f ): 



dX? 
dt 



—tt — ^ K ab (Xi — Xi) 



b 



with 



K 



m b AD a D b \VaW ab 



(17) 



(18) 



PaPbD a + Db \r ab \ 

where the subindexes a and b are used to denote different 
SPH particles. As a first approximation, in this paper we 
have considered a constant diffusion parameter D a — D b — 
D for all the gas particles. Except in the first three synthetic 
tests of § 14.11 where the instantaneous mixing of metals 
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10 7 
log 10 (t) (yr) 

Figure 1. Composition of the ejecta as a function of time for three SSPs with solar total metallicity but different abundance ratios. The 
solid line corresponds to solar proportions, the dashed line to [O/Fc]=0.5, and the dashed-dotted line to [O/Fe]=-0.5. It can be seen how, 
even at identical total metallicity, in many cases the production of a particular isotope varies significantly. This effect is due solely to the 
use of the Qjj matrix formalism, since the stellar evolution models used in this work only take into account the initial total metallicity, 
which is the same for the three compositions. 



through the gas component is assumed, and in the diffusion 
test of § 14.1.41 we have adopted the value D = 9.25 x 10 26 
cm 2 /s as suggested by LPC02. 

In our code, the diffusion equation (|17p is computed 
at each timestep by considering all the active gas parti- 
cles a, and their active neighbours By active particles 
at timestep [t,t + At], we mean those needing the update 
of their physical properties at that timestep. In order to 
solve Eq. (|17|l it must be noted that metal diffusion would 
introduce a new time-scale (X/X) jjg- and therefore, an ad- 
ditional criterion limi ting the time stepp ing. Just like for the 
radiative cooling (see lSerna et al,l l2003), such an additional 
At-control criterion can be circunvented by solving Eq. (|17|) 
in integrated form. We have then used the fact that, because 



of the Courant condition, the density field (and hence K a ) 
is nearly constant over a timestep. Therefore, Eq. (|17p can 
be anallytically integrated to find: 



AX? = X? (t + At) - X? (t) = (W) 

b 



with 



AXf = -AX? a = i (l - e - 2KabAt ^j (Xt(t)-X?(t)) .(20) 



4 Only active gas particles are considered in order to maintain the 
symmetry and ensure the abundance conservation in the diffusion 
process 



Such analytical expressions are those actually used in our 
implementation of the diffusion process. 
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3 THE COMPOSITION-DEPENDENT 
COOLING FUNCTION 

3.1 The cooling model 

The publicly available MAPPINGS III code (see 
ISutherland fc Dopital Il993l , hereafter SD93, for a much 
wider description) consists of a detailed cooling model for 
low and high temperatures, which includes calculations for 
up to 16 atoms (H, He, C, N, O, Ne, Na, Mg, Al, Si, S, CI, 
Ar, Ca, Fe and Ni) and all stages of ionization. In this code, 
the net cooling function of an optically-thin astrophysical 
plasma is obtained by adding the contribution of different 
processes: collisional line radiation (including fine-structure, 
inter-system, and forbidden emission), Aij ncs , free-free and 
two-photon continuum radiation, A cont , recombination 
processes with both cooling and heating effects, A rC c, 
photoionization heating, Aphoto, and collisional ionization, 

Acoll- 

A ne t Ai; nes ~t~ A C ont 

± A roc — Aphoto + Acoll (21) 

All the detailed cooling computations in this paper have 
been carried out by using MAPPINGS III. To obtain results 
easily comparable to those reported by SD93, widely used 
in cosmological simulations, we selected the algorithms and 
assumptions of its previous version. In addition, all these 
computations considered Collisional Ionization Equilibrium 
(CIE) conditions. Therefore, the photoionization heating 
is insignificant, leaving only collisional line radiation, con- 
tinuum radiation, and recombination heating as significant 
terms. 

3.2 Metallicity dependent cooling in cosmological 
simulations 

In principle, the accurate value of the normalised cooling 
functioqj depends on both the local temperature and the 
detailed chemical composition: 

Ajv=Ajv(T,Z) (22) 
where 

Z = {Z 1 ,Z 2 ,...,Z N ) (23) 

is a vector containing the abundance Zi of the N = 16 chem- 
ical species taken into account by MAPPINGS for the cool- 
ing rate. In a cosmological simulation, one must deal with a 
huge number of gas particles with different metal mixtures. 
A full computation (e.g., by directly using MAPPINGS) of 
the cooling function of each particle according to its detailed 
chemical composition would be very expensive in terms of 
computing time. 

To deal with the above difficulty, the most common ap- 
proach consists of reducing the Z dependence of the cooling 
function into a much smaller number of parameters or, ide- 
ally, into just one parameter f. 

Aj V =Aj V (T,C) (24) 

5 Please note that we define the normalised cooling rate as 
Ajv = A ne ^/p 2 , instead of Ajv = A ne ^/n e nt, where n e is the 
electron number density and nt the total ion number density. 
The resulting units for Ajy are erg cm 3 s _1 g — 2 . This approach 
eases the computation of A ne ^ in SPH codes. 




^ 0.02 0.04 0.06 0.08 0.1 



Figure 2. Cooling function dependence on the total metallicity 
Z for a given (T = 10 5 5 K) temperature for distribution of abun- 
dances present on a cosmological self-consistent simulation with 
metal enrichment. The solid line represents the cooling according 
to SD93. 

For example, SD93 have considered the total metallic- 
ity, Ztot, as their choice of f. More specifically, the cooling 
function is computed and tabulated by using chemical pro- 
portions interpolated, according to the total metallicity Ztot , 
between primordial and solar abundances. Such an approach 
allows for a fast implementation of a metallicity dependent 
cooling function in cosmological simulations. 

In order to test the above approach, we have considered 
a sample of gas particles with different metal mixtures. In- 
stead of assigning a random composition to each particle, we 
have randomly selected ~ 10 particles extracted from cos- 
mological ACDM simulations with the chemical evolution 
model described in § [2] 

A full computation of the cooling function has been 
carried out for each particle by using MAPPINGS and its 
individual chemical composition. The dots in Fig. [5] repre- 
sent the individual values of the cooling function when the 
same temperature (T = 10 5 ' 5 ) is assigned to all the gas par- 
ticles in the sample. The solid line represents instead the 
cooling rate obtained from the SD93 method. It can be seen 
from this figure that the latter approach gives a reasonable 
approximation of the cooling rate at different Z to t values. 
However, the figure also shows an important dispersion on 
the cooling rate of gas particles with the same total metal- 
licity but different metal mixtures. Such a dispersion could 
lead to errors in the estimate of At(Z) of almost one order 
of magnitude for sub-solar metallicities. 

Different approaches can be envisaged to improve the 
modelling of the cooling function in cosmological simula- 
tions. For example by characterising the composition de- 
pendence of the cooling function with more than one pa- 
rameter. This could be achieved by considering, in addition 
to the total metallicity, the alpha-element enhancement or 
any other parameter providing a more detailed description 
of the metal content. Another approach that needs to be 
explored consists of maintaining the Eq. (|24l) . with just one 
metallicity parameter, but using a different choice of £ to de- 
scribe the effect on the cooling function of different chemical 
compositions. 
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Figure 3. Cooling function dependence on the fx parameter for 
a given (T = 10 K) temperature. 



In order to analyse this latter possibility, we have em- 
ployed a Dimens ion Reduction Regression (DDR) technique 
(|Weisberd I2002T ). In such a procedure, one tries to reduce 
the multidimensional dependence of a function A(Z) (e.g., 
the cooling function at a given temperature) to a small 
number d of parameters expressed as linear combinations 
fi = ci ■ Z, C,d = Cd ■ Z, where c, are vectors, of the 
same dimension as Z, containing the coefficients of each lin- 
ear combination fj. If d is very small, one or two, then the 
regression problem can be summarised using simple graph- 
ics (a single A(£i) plot for d = 1, or a 3D plot for d = 2) 
containing all the regression information. 

Several methods for estimating d and the relevant coeffi- 
cients ci, Cd have been suggested in the literature. In this 
pape r we have used a sliced inverse regression method (0 
where the range of At values is divided into several 
intervals, or slices, and then a weighted principal component 
analysis is performed. This gives higher importance to the 
slices with higher covariance. The eigenvectors of the covari- 
ance matrix, ordered by their corresponding eigenvalues, are 
the preferred directions (i.e. the d Cj vectors). 

When the above algorithm is applied to the cooling 
function, one finds that just one parameter 



C(T) = c(T) • Z 



(25) 



is enough to accurately fit Am at a given temperature T. For 
the sake of clarity, Table 1 gives the resulting c(T) coeffi- 
cients for some temperature values a nd for the same metal s 
as those considered in the yields of iGavilan et al. | (|2005h . 
The metallicity parameter in Table 1 has been normalised 
according to ( N (T) = [C(T) - Co(T)]Ci(T) so that ( N (T) 
lies in the interval [0, 1]. The corresponding (o(T) and Ci(^) 
normalisation parameters are also given in Table 1. An ex- 
tended version of this Table is available online. 

To check the above result, we have repeated the same 
test as in Fig. [5] For each gas particle in the sample we have 
considered its detailed chemical content to determine its cor- 
responding ( parameter (Eq. I25|l , as well as to carry out a 
full MAPPINGS III computation of its individual cooling 
rate at T = 10 5 ' 5 K. The dots in Fig. [3] show the {Aj,(j) 
values obtained for each particle j. It can be seen that such 




T(K) 

Figure 4. Cooling function for a sample of particles with so- 
lar metallicity (Z = 0.02 ± 0.0002) but different metal mixtures 
as found in a cosmological simulation (thin lines). The thick 
black line represents the cooling function for solar abundances 
i Anders fc Grevessdl 19891) . 



a plot now presents almost no dispersion and, therefore, the 
Ajv(T, £) values can be easily and accurately tabulated (see 
Table 2). 

Summarising, given a gas particle with known temper- 
ature T and chemical content Z, its appropriate metallicity 
parameter ((T, Z) can be computed by using Eq. (I25|l and 
the coefficients of Table 1. An accurate estimate of the cool- 
ing rate can be then obtained through interpolation from 
Table 2. The resulting algorithm needs an almost negligible 
amount of computing time but implies a remarkable im- 
provement on the cooling modelling, as compared to other 
approaches based on Z to t ■ Indeed, we have applied this algo- 
rithm to a sample of simulated gas particles with Z to t ~ Zq 
but different chemical contents. Fig. [4] shows their corre- 
sponding cooling functions (thin lines), as well as the cool- 
ing function (thick line) obtained for solar abundances (i.e., 
in the case considered by SD93). 

The physical meaning of the c(T) values obtained in 
the above DDR procedure can be understood as coefficients 
giving, at a given temperature T, the weight of each element 
on the total cooling function. For example, Fig. [5]-a shows 
the logarithmic partial derivative of Ajv with respect to the 
abundance of each element, i.e., a quantity roughly giving 
the weight of each element on the total cooling function. 
Fig. [5]-b shows instead the contribution = CiZi of each 
element i to the metallicity parameter £(T). In both panels, 
solar proportions have been assumed. It can be seen from 
this figure that both quantities are correlated, with the most 
important coolants contributing more to the metallicity pa- 
rameter C(T). 

It is important to note that, in SPH simulations, specific 
internal energy u is tracked instead of temperature T. Both 
are related by 



Ut(Z) = 



2fim p 



-k B T . 



(26) 



is the proton mass 



where ks is Boltzmann's constant, 
and fj, is the mean molecular mass. For simplicity, some SPH 
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T(K) 

Figure 5. (a): portion of the total cooling function contributed 
by C (solid line), O (dashed line), Ne (dashed-dotted line), Mg 
(dotted line) and Fe (dashed-triple dotted line) for a solar mixture 
of elements (Z = Zq). (b): corresponding values of the CiZi 
terms of C,t- 



codes fix fi to a constant value regardless of metallicity. In 
some cases (e.g., iKawata fe Gibson! 120031 ) such a constant 
value is chosen to represe nt a fully ion ised gas (ft = 0.6), 
while in other cases (e.g.. iBerczikl fl999h the adopted value 
represents a cool gas with primordial (ft — 1.2) or solar 
(ft = 1.3) abundances. To be consistent with our aim of de- 
veloping a model that takes into account the full dependence 
on the chemical composition, we have preferred to consider 
the mean molecular mass as a function of T and Z. Conse- 
quently, we have applied a DDR procedure similar to that of 
Eq. (|24|l to write u(T, Z) as a linear combination of Z. The 
corresponding DDR coefficients are also available online. 



4 TESTS AND RESULTS 
4.1 Synthetic tests 

In order to test the results of our model, as well as the res- 
olution effects, we have carried out a series of tests based 
on classical evolution models for chemical evolution. In all 
these tests, particles are just discrete mass elements submit- 
ted to certain constraints. Therefore the models presented in 
this section are useful to test the different aspects of our im- 



plementation of the metal enrichment, uncoupled from any 
dynamical effect. 

4-1.1 A single stellar population 

Following LPC02, a first basic test consists of analysing the 
chemical evolution in the simple case of a single burst of 
star formation. Such a test is designed to check the validity 
of our statistical approach by comparing its results to the 
expected chemical production of an SSP. 

The single-burst model begins with N = 5000 newly 
born star particles of total mass 10 11 M©. All particles have 
primordial composition at t — and are left to evolve after- 
wards according to the statistical prescriptions for gas resti- 
tution and metal enrichment described in §[2] Consequently, 
as time progresses, some stellar particles will turn into gas. 
At each time step, the metal content of the gas component 
is enriched through the yields of the remaining stellar parti- 
cles according to either Eq. or (I12|l . In this test, metals 
are instantaneously mixed through the gas component and 
no further SF episodes take place. 

Fig. [6] displays, as a function of time, the gas mass frac- 
tion and the chemical composition of the gas. The solid 
lines correspond to the exact analytical predictions for an 
SSP (i.e., those obtained in the continuous limit by directly 
integrating Eqs. 151 and [T2l for a unique SSP). The dashed 
and dotted-dashed lines give instead the results obtained 
from a statistical model of chemical evolution based on Eq. 
(|12|) and Eq. respectively. It can be seen from this 

figure that both statistical methods closely reproduce the 
expected trends for an SSP. Therefore, our scheme of pro- 
gressive metal ejection (Eq. I12|l gives, within a stochastic 
model, a fair representation of the overall production of an 
SSP. 



4-1.2 Closed box model of an elliptical galaxy 

Another useful te st for chem ical evolution models is the 
closed box model (lPagellll997l ) with one zone. We have im- 
plemented such a test by initially considering N gas parti- 
cles that, during the simulation, randomly turn into stars 
according to a pre-defined probability. More specifically, in 
order to mimic an elliptical galaxy, we have considered an 
exponentially decaying star formation rate: 



(27) 



so that, assuming no gas feedback, the rate of change of the 
gas mass fraction, g(t), is given by 



dg(t)/dt = . 



(28) 



The above equation implies that the probability p that 
a gas particle forms stars in a time At is 

f„—T~At 



1 K 



p = 



(29) 



K + e rt (T - k) 

where k can be written in terms of the final fraction of gas, 
9i- 



■ T (gi - 1) 

1 - e T 



(30) 



As in the previous test, all particles are left to evolve 
according to the statistical model for gas restitution and 
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metal enrichment described in § [2j except that here we use 
Eq. (|29[) to compute p at each time step for all gas par- 
ticles and draw random numbers to decide which particles 
actually form stars. We again consider that metals are in- 
stantaneously mixed through the gas component. 

The above closed box model then gives a simple repre- 
sentation of objects with a more complex sequence of star 
formation bursts and, therefore, containing a mixture of dif- 
ferent SSPs. The main advantage of this model is that it has 
a known theoretical prediction for the resulting metal distri- 
bution function (MDF) of the gas component (see Appendix 
[5]). The comparison with such a theoretical MDF then con- 
stitutes a strong test, uncoupled from hydrodynamic and 
gravitational effects, of the stability of particle-based nu- 
merical methods against a degraded resolution. 

We have run this test for four different particle numbers 
(N = 10000, 1000, 200 and 100). In all runs we considered 
t — 1/5 Gyr - , a typical value for massive elliptical galaxies, 
and a constant timestep of At = 1.38 10 6 yr (i.e., 10000 in- 
tegration steps over a Hubble time). In order to have enough 
resolution for the gas component, we have used a high value 
(gi — 0.4) for the final gas fraction. In a series of runs we 
have used Eq. f) to incorporate the stellar production of 
metals into the gas component, whereas in another series of 
runs we have used the new procedure proposed in this paper 
(Eq. I12p to account for the metal feedback. Fig. [7]shows the 
results of this test, where the solid line corresponds to the 
theoretical MDF, whereas the particle-based predictions are 
shown as dashed (for Eq . 112(1 and dotted-dashed (for Eqllip 
lines. As can be seen from Fig. [7] for a high number of parti- 
cles (N = 10000) both procedures for metal feedback lead to 
MDFs in good agreement with the theoretical solution. Ob- 
viously, the kernel procedure used to draw the particle-based 
results leads to MDFs where the sharp peak and cutoff of the 
theoretical solution at [O/H] ~ 0.23 are smoothed. In tests 
with N = 1000 particles, both metal feedback procedures 
still give results in reasonable agreement with the theoreti- 
cal solution, although the statistical noise is slightly less vis- 
ible in the predictions obtained from Eq. (|12[1 than in those 
found from Eq. (|11|) . At low particle numbers (N — 200 and 



100), the results from Eq. (|12[1 degrade sensibly better and 
keep the shape of the MDF almost unaltered, while those 
from Eq. pip become noise-dominated due to the lack of 
enough sampling data. 



4-1.3 Multi-zone model of a Spiral Galaxy 
In order to test our implementation of the Qg formalism. 



we ha ve also performed a test based on the IMolla fc Dfaj 
(2005) multi -zone evolution model for a disk galaxy. Such a 
model also considered the Qij formalism and used the same 
stellar libraries as in our code. To e nsure that our results 
are directly comp a rable to those of IMolla fc Dfaj (|2005l ) 
and Gavilan et all (12005). we have adopted the same IMF 
(IFerrini et al. I ll990l , ll992h for this test as that used in such 
studies. 

The lMolla fc Di'aj (|2005l ) model considers a galaxy with 
a spherical halo of radius Rh and a concentric cylinder con- 
stituting the disk. The spherical halo is divided into concen- 
tric cylindrical regions 1 kpc wide with a height determined 
by the corresponding galactocentric distance on the disk and 
the total radius of the sphere. The corresponding regions on 
the disk are also concentric cylindrical shells 1 kpc wide, but 
with a constant height ho ~ 0.2 kpc. 

In the initial conditions, the halo contains the total mass 
of the galaxy in gas phase. The halo mass has a radial dis- 
tribut ion consistent with the rotation curve of I Persic et al.l 
(1996). In units of 10 9 Mq, such a distribution is given by 
M(R) = 2.32 x 10 s RV{Rf, with 



V(R) = V opt 



-, 1/2 



0.72 



(x 2 +0.61) 1 - 



+ 1.07 



x 2 + 2.25 



.(31) 



Here, x — R/R opt and R op t = -Rsr/2.5. We have con- 
sidered V opt = 200 km/s, R op t = 13 kpc and L = W 10 4 Lq, 
which correspond to galaxies similar to the Milky Way. 

For each region i the matter in the halo can be ei- 
ther in the form of stars of diffuse gas, with total mass 
M s h and M g h, respectively. In the corresponding disk re- 
gion, the model allows for the following phases: diffuse gas 
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Figure 7. Metallicity Distribution Function (MDF) for the stellar particles in the closed box test. The different panels correspond to a 
decreasing number of particles (TV =10000, 1000, 200 and 100) and, therefore, to a decreasing mass resolution. In order to make easier 
the comparison of the different curves, the particle-based MDFs are displayed as ker nel-smoothed distributions with the bandwidth, 
h = 1.06aN~ 1 ' 5 , where a is the standard deviation of the sample, as suggested by IScottl l| 19921) . The solid line corresponds to the 
analytical solution, the dashed-line to the numerical results obtained from Eq. JT2}, and the dotted-dashed line to those found from Eq. 

CC0- 



(M g d), clouds (M c d), and stars (M s( j). In this latter phase, 
the adopted IMF implies a certain mass M S 2 of stars with 
M > 8 Mgfl The mass of the different phases can change 
through the following conversion processes: i) Diffuse gas in- 
fall from the halo to the disk, with a rate rhiMF, ii) In the 
halo, star formation from diffuse gas (itisfh), hi) In the disk, 
star formation from clouds (rhsFu), either from cloud-cloud 
collisions and from massive star-cloud interactions, and iv) 
Cloud formation from diffuse gas in the disk (tticfd)- The 
complete set of equations for each zone is therefore (see 
iMolla et al.1ll99fj , for details): 

dMg h /dt = W H - rhiNF - itisfh 
dM sh /dt = rhsFH ~ W H 

6 Note that the mass of massive stars in an SSP, M s 2, declines 
with the SSP age and eventually drops to zero. 



dM sd /dt = m S FD-W D (32) 
dM cd /dt = mcFD-rhsFD 
dMgd/dt = rhiNF - tticfd + Wd 

where Wh and Wd represent the gas return rate for the halo 
and disk, respectively, and 

miNF = M gh /r 
msFH = KMgj? 

mcFD = (33) 
rhsFD = HMc d + e a M cd M s2 

with r = r c exp[{R - R c )/Xn], K = e K (G/V H ) 1/2 , V = 
e^G/Vo) 1 ^, and H = e„/V D . Here, R c = R opt /2, X D = 
0.15R op t, G is the universal gravitational constant, Vh is 
the volume of the halo region and Vd is the volume of the 
disk region. We have considered parameter values that cor- 
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respond to a galaxy mimicking the Milky Way: r c = 4 Gyr, 
e K = 5.3 10" 3 , e M = 0.80, e H = 0.28 and e a = 0.83. 

The above model has been implemented by considering 
16 regions in the disk, from R = 2 kpc to R — 18 kpc. We 
explicitely avoid the bulge since our model is not applicable 
there. Each region is sampled with 



N = 1000 max 



1. 



1 



1 — exp( — 1/t) 



(34) 



particles of individual mass M(R)/N. All particles are ini- 
tially labelled as diffuse gas in the halo so that, except 
for M g h, all phases have a vanishing total mass. At each 
timestep At, the rate of change for the mass of the different 
phases is computed using Eqs. (|32[) -(133 [) . The corresponding 
change on the mass Mi of each phase i is then computed 
at the first level of approximation, AMi = (dMi/dt)At, and 
then expressed as discrete changes in the numbers of parti- 
cles. After re-labelling the particles, our routines for chem- 
ical evolution are called to compute the metal production 
and return of gas (terms Wh and Wd in Eqs. I32|l in both 
the halo and the disk, except for the star formation, which 
is computed in this text from Eq. (|33l) . 

Fig. [5] c o mpar es our results and those found by 
iMolla fe Di'az] (|2005l ) at t = 13.72 Gyr. We see from this 
figure that the radial distribution of the different elements 
and model components o btained from o u r code are in close 
agreement with those of IMolla fe Dfa3 l|2005h . The small 
deviation of the oxygen and carbon abundance at the outer 
galactic regions is probably due to the fact that our approach 
is based on masses that are sampled by a discrete number 
of particles. Both methods give also very similar results for 
the time evolution of the oxygen and carbon abundance in 
the gas for the solar cylinder. 

Fig. [9] shows the time evolution of the stellar [O /Fe] 
ratio obtained from the Qij formalism (dashed line) as com- 
pared to that found when the Qij formalism is not used 
(dashed-dotted line). This latter result has been computed, 
for each stellar particle with composition Xj, by using in 
Eq. @ effective abundances X'j scaled, assuming solar pro- 
portions, according to its total metallicity. We see from this 
figure that the latter procedure underestimates the [O/Fe] 
ratio. This is due to the fact that forcing solar proportions 
for a stellar particle with [O /Fe] > implies the use in Eq. ((9} 
of a larger effective abundance of iron (X Fc > Ap c ). Con- 
sequently, when computing the ejection of iron eFe(i), the 
dominant affected term in Eq. ((9} is QFe,FeA" Fe > ^Fe.Fe^Fe, 
and eFc(t) is overestimated. Although less significantly, sim- 
ilar arguments imply that the ejection of oxygen is underes- 
timated. These two effects can also be noticed in Fig. [T] and 
result in the underestimation of the [O/Fe] ratio observed 
in Fig. E 



usion in a homogeneus box 



4.1.4 



The diffusion of metals in our model remains to be tested. 
In order to test our implementation of Eq. (|16p . described 
in § 12.4.21 we have devised a simple test. Starting with an 
isolated cubic box of L — 1.43 Mpc filled with 32 3 particles, 
distributed on a regular grid and representing a cloud of 
homogeneous gas with density p = 1.058 x 10~ 35 g/cm 3 
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Fi gure 8. Comparison of our results (thin lines) and those found 
by IMolla fc Diaa ||2005| ) (thick lines) in their multi-zone model 
of a spiral galaxy, (a) Top panel: radial distributions of atomic 
gas, molecular clouds and stellar density at the end of the run. 
Bottom panel: radial distributions of oxygen and carbon, (b): 
Time evolution of oxygen and carbon abundances in the gas for 
the solar cylinder. 




Time(Gyr) 

Figure 9. Time evolution of the [O/Fc] ratio of stellar parti- 
cles in the solar cylinder. The dashed line corresponds to the 
result obtained from the Qij formalism. The result displayed as 
a dashed-dotted line has been obtained by computing the metal 
production of each stellar particle using its total metallicity and 
assuming solar proportions. 
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Figure 10. Upper panel: final spherical distribution of metallicity 
for the diffusion test (crosses) and analytical solution (solid line) 
for t = 13.72 Gyr. Lower panel: time evolution of the metallicity 
at r = 200 kpc. 



and let diffusion algorithm act for 13.72 Gyr. A diffusion 
constant of D = 4.63 x 10 29 cm 2 /s has been used, and no 
other processes have been considered, so that the particles 
are kept fixed to the grid. 

This simple setup allows for an analytical solution of 
the diffusion equation with initial conditions 

X(r,t = 0) = 0.LS(r) , (35) 

namely 

X(r,t) = / J^^^^e"^*^- . (36) 

While initial conditions of Eq. (135[) cannot be imposed in a 
discrete model, both should converge if r or t are big enough. 

In fact, the upper panel of Fig. HOI shows the comparison 
between the final metallicity profile of the numerical test and 
its corresponding analytical solution. We find a good agree- 
ment for such a high t. Lower panel shows the metallicity 
time evolution in a bin centered around r = 200 kpc. We 
find the numerical solution to have a shallower growth than 
the analytical one, due to the fact that the numerical initial 
conditions are not an exact delta function. At later times, 
however, both solutions converge. 



4.2 Cosmological simulations 

We have finally performed different simulations of galaxy 
formation and chemical enrichment within a cosmological 
context. All simulations started at redshift z = 20 from 
initial conditions that are a Montecarlo realisation of the 
field of primordial fluctuations to a concordance cosmolog- 
ical model (a flat ACDM model, with h = 0.7, £Ia = 0.7, 
Qm = 0.3 and fib = 0.04). The as value has been taken 
slightly high (<t $ = 1) in order to mimick an active region 
of the universe l|Evrard et alJll990h . The evolution of these 
fluctuations was numerically followed up to z = by means 
of th e parallel version of the DEVA code (see ISerna et al.l 
2003, for a detailed description of DEVA) where we have 
coupled our model for chemical evolution and cooling. In the 
test performed in this section we have used a star-formation 
threshold of p t h — 3 x 10 -25 g/cm 3 , an efficiency of c» = 0.3, 
a gravitational softening of e = 1.5 kpc//i and N = 2 x 64 3 
particles in a box of L = 7 Mpc/fr, implying a mass res- 
olution of nib = 2.06 x 10 7 M© for baryonic particles and 
m-dm = 1-34 x 10 8 Mq for dark matter particles. 

Individual galaxy-like objects of different morphologies 
naturally appear in these simulations as a consequence of the 
cosmic mass assembly and the physical processes taken into 
account. Differently from the semi-analytical models of § 14.11 
no assumptions are made about the spatial variations of the 
infall rate of gas and its relative importance as compared to 
the star formation rate. 

These simulations then constitute an appropriate test 
to analyse: i) whether the resulting galaxy-like objects have 
realistic chemical properties; ii) the possible effects of in- 
cluding the cooling function presented in § [21 that takes 
into account the full dependence on the detailed chemical 
composition of gas particles; and iii) the importance of us- 
ing the Qij formalism in cosmological simulations. In or- 
der to address these issues we have carried out different 
runs. In a first run (referred as the ^"-cooling run) we used 
the composition-dependent cooling of § [3] whereas in a sec- 
ond run (referred as the Ztot-cooling run) we used the total 
metallicity-dependent cooling of SD93. In both cases, we 
used the Qij formalism of § [5] In a third run (referred as 
the no-Qij run), the metal production of each stellar parti- 
cle was instead computed by using its total metallicity and 
assuming solar proportions. 

Elliptical-like objects (ELOs) const itute a family of ob- 
jects with simple scaling relations (see lOnorbe et al. I l2005l . 
for a study of the origin of the Fundamental Plane of ELOs 
obtained in DEVA simulations). In this work, we focus on 
the study of ELOs because given their simplicity, less reso- 
lution is needed to properly simulate them. ELOs were iden- 
tified as those objects having a prominent stellar spheroidal 
component with hardly disks at all. Here we will focus on 
the analysis of the chemical properties of the most massive 
ELOs at z — 0, sampled with at least 1000 baryonic parti- 
cles. This selection criterion produces eight massive ELOs, 
with a stellar mass range of ~ 3 x 10 10 - 1.5 x 10 11 M Q . 

In both the ^"-cooling and Z to t-cooling runs, we found 
that ELOs have metallicities with individual mean values 
that are consistent with t hose observed for ell iptical galax- 
ies with similar masses l|Thomas et all l2005h . For exam- 
ple, Fig. [TT] shows the age-metallicity relation and metallic- 
ity distribution function (MDF) for the most massive ELO 



14 F.J. Martinez- Serrano, A. Serna, R. Dommguez-Tenreiro and M. Molld 




[Fe/H] 

Figure 11. Chemical distribution for the most massive ELO in 
our simulation. The upper panel shows the age-metallicity rela- 
tion as a contour plot. Only the (-cooling run is displayed since 
there are virtually no differences between the relations arising 
from both cooling methods. The lower panel shows the stellar 
metallicity distribution function (number fraction of stars with 
a given metallicity) obtained in the (-cooling (dashed line) and 
^tot -cooling runs (das hed-dotted line) for RGB stars in the halo 
(seeH 

arris et ah 2007 for an observational counterpart). 




M»(M ) 

Figure 12. Inverse concentration index C~ 1 = r p 5o/r p go as a 
function of total stellar mass of the ELOs found in the Ztot- 
cooling (circles) and (-cooling (stars) runs. 



dius to the 90% mass Petrosian radius). It can be seen from 
this figure that the Z to t-cooling method produces objects 
that are more concentrated and with a higher stellar con- 
tent than those found in the ("-cooling run. Such differences 
are small for the most massive objects but become impor- 
tant as we consider less massive ELOs. They are probably 
due to the fact that the Z to t method has some tendency to 
overestimate the cooling rate of particles not having solar 
abundance ratios (see Fig. [2}. 

The effects of using the Qij formalism in cosmological 
simulations are shown in Fig. 1131 Such a figure displays, for 
the most massive ELO, the time evolution of the [O/Fe] ratio 
of stellar particles with radial distances r < 15 kpc. The 
dotted-dashed line corresponds to the result obtained from 
the no-Qij run, whereas the dashed line displays the result 
obtained when using the Qij formalism and the same cooling 
method as in the no-Qij run (i.e., the (-cooling run discussed 
above). We see from this figure that, just like in the multi- 
zone test of § 14.1.31 the assumption of solar proportions leads 
to a significant underestimation of the [O/Fe] ratio. 



(M = 1.5 x 10 M ) identified in such simulations. We have 
found that the age-metallicity relation obtained for this ob- 
ject in the (-cooling (shown in the upper panel of Fig. Ill[) is 
almost indistinguishable from that found in the Z to t-cooling 
run. Such a relation shows a very fast enrichment on the first 
3 Gyr, when most of the star formation happens, followed 
by a very slow enrichment at recent times. The large scatter 
at t < 3 Gyr is probably due to the fact that stars were 
formed on separate smaller objects that merged later on to 
form the final object. The MDFs found for this object in 
both simulations (lower panel of Fig. Ill|l are also very sim- 
ilar and with a shape that clo sely resembles that observed 
for some giant ellipticals (e.g., Harris et al.ll2007l ). In addi- 
tion, the central oxygen abundance (not shown in Fig. Ill[) of 
this ELO has a value ([0/H]~ +0.2) that agrees with that 
expec ted from the mass-me tallicity relation given by the eq. 
(3) of lThomas et alj i|2005l ) 

In less massive ELOs, some significant differences ap- 
pear between the results obtained from the ( and Ztot 
cooling runs. Fig. 1121 shows the inverse concentration index 
C _1 = r P 5o/r P 9o (i.e., the ratio of the half-mass Petrosian ra- 



5 SUMMARY AND CONCLUSIONS 

In this paper we have introduced an SPH model for chemical 
enrichment and radiative cooling in cosmological simulations 
of structure formation. Particular attention has been paid to 
including, by means of fast algorithms, the full dependence 
of both processes on the detailed chemical composition of 
star and gas particles. 

As compared with previous implementations in N-body 
simulations codes, the main features of our SPH-model for 
chemical evolution and cooling are: 

i) Our model takes into account the delayed gas resti- 
tution from stars by means of a probabilistic approach that 
shares many aspects with that proposed by LPC02, except 
that the stellar yields are progressively distributed through 
the neighbouring gas (Eq.ll2[). The tests described in § 14.1.11 
and 14.1.21 show that our scheme, as well as that of LPC02 
(based on Eq. Hip , provides a fair representation of the 
chemical production of an SSP and reproduces the expected 
trends for both one-zone and multi-zone models of galax- 
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Time (Gyr) 

Figure 13. Time evolution of the stellar [O/Fe] ratio in the most 
massive ELO. The dashed line corresponds to the result obtained 
from the Qij formalism, whereas the dotted-dashed line displays 
the result obtained in the no-Qij run. In both cases, the (-cooling 
method was used. 



ies. Nevertheless, a scheme based on our Eq. (|12|l reduces 
the statistical noise and, at low particle numbers, degrades 
sensibly better than one based on Eq. (|lip , which becomes 
noise-dominated due to the lack of enough sampling data. 

ii) The full dependence of the metal production 
on the detailed chemical composition of stellar particles 
is taken into account by means of the Qij formalism 
(Talb o"t~fe Arnettl |l97ot ). that relates each nucleosynthetic 
product to its sources. Therefore, the assumption of solar 
proportions is relaxed in our model. Although the Q^- for- 
malism has b een previously considered in chemical ev olution 
models fe.g.. lFerrini et al.lfl99^ ; iGavilan et al.ll2005t) . this is 
the first time that it is used in N-body simulation codes. Re- 
leasing the assumption of solar proportions is important to 
accurately follow the evolution of abundance ratios like, e.g., 
the enhanced [a/Fe] ratio observed in elliptical galaxies and 
in spiral galaxy bulges. This ratio, as well as its dependence 
on the central velocity dispersion, could provide us with im- 
portant constraints on galaxy formation models. Indeed, the 
tests performed in this paper for both multi-zone models (§ 
I4.1.3|l and cosmological simulations f§ I4.2|l suggest that the 
assumption of solar proportions leads to a significant un- 
derestimation of the [a/Fe] ratio in simulated galaxy-like 
objects. 

iii) In the same way, the full dependence of radiative 
cooling on the detailed chemical composition of gas parti- 
cles has been implemented through a fast algorithm based on 
a metallicity parameter, ((T), that takes into account the 
weight of the different elements on the total cooling func- 
tion. We have compared the results obtained when such a 
composition-dependent cooling is used (referred as (-cooling 
runs) and those found from the total metallicity cooling ta- 
bles of SD93 (referred as Z to t-cooling runs). To that end, 
we have carried out (§ 14. 2[) different simulations of galaxy 
formation in the framework of a concordance cosmological 
model. For the most massive elliptical-like objects (ELOs), 
we have found that the differences between the ( and Z to t- 
cooling runs were small. Such massive ELOs are charac- 
terised by strong and short bursty star formation events 
at early times, where most gas is exhausted and, therefore, 



where dynamics soon becomes dissipationless. However, for 
less massive ELOs, some important differences appear be- 
tween the results obtained from the ( and Ztot-cooling runs. 
Probably due to the fact that the Z to t method has some ten- 
dency to overestimate the cooling rate of particles not having 
solar abundance ratios, such a method produces ELOs that 
are more concentrated and with a higher stellar content than 
those obtained in the ("-cooling run. Such differences become 
larger as we consider ELOs with smaller stellar masses. 

The above scheme for chemical enrichment and cool- 
ing has been i mplemented in t he parallel-OpenMP ver- 
sion of DEVA |S crna et al.ll2003l ). a Lagrangian code par- 
ticularly designed to study galaxy assembly in a cosmo- 
logical context. In this code, gravity is computed through 
an AP3M-like method, while hydrodynamics is computed 
through an SPH technique with algorithms and correction 
terms ensuring an accurate implementation of conserva- 
tion laws (energy, entropy and angular momentum). The 
DEVA code has already been abl e to produce both elliptical 
ilDommguez-Tenreiro et alj|2004 l200rj; lOfiorbe et al.ll2005l 
l200rj . |2007l ) and spiral (|Saiz et alj|2002l ) galaxy-like objects 
with realistic structural, kinematic and dynamical proper- 
ties. Using this model for chemical evolution and cooling, we 
will analyse in forthcoming papers some important chemi- 
cal properties like, e.g., the metallicity distribution func- 
tions (MDF) and abundance gradients within individual ob- 
jects of different morphologies and environments, as well as 
the fundamental metallicity relations (mass-metallicity and 
mass- [a/Fe] ratio) for large samples of ELOs obtained in 
cosmological simulations. 
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APPENDIX A: ANALYTICAL MDF FOR A 
CLOSED BOX MODEL WITHOUT 
INSTANTANEOUS RECYCLING 



Following iPagell (|l997f ). we consider a closed box filled with 
primordial gas. We assume a IChabrierl (|2003ah IMF, 4>(M), 
and the star- forming law given by Eq. (|27[1 . We denote 
the gas and star fractions by g(t) and s(t), respectively. The 
mass conservation then implies that 



g(t) + s(t) = 1 



and therefore 
ds dg 
crt ~ di 
where e(t) is the ejection rate of gas: 



*(*) " <t) 



e(t) = 



(M-M r (M))*(i- 



M T (t) 



T (M))^P-dM 

y " M 



(A-l) 



(A-2) 



(A-3) 



Here t(M) and M T (t) represent the lifetime of a star of 
mass M and its inverse, respectively. 

Finding the MDF is equivalent to expressing the star 
formation rate as a function of the metallicity. We must 
then obtain the time dependence of the gas composition, 
Z(t), and invert such a function to find t(Z). By substituting 
the latter function into the star formation rate Vt(i), we will 
have *(Z). 

To find Z(t), we notice that the time derivative of the 
total content of metals g(t)Z(t) must be equal to the total 
ejections by dying stars ez(t) minus the mass of metals that 
gets trapped in newly formed stars "if(t)Z(t): 

d(Zg) 



dt 



e z (t) - *(t)Z(t) 



(A-4) 

Using Eq. IA-21 we can simplify this expression: 

e(t)Z[t) . (A-5) 



9(t)^=ez(t) 



The expression for ez(t) is similar to that of e(t): 
ez(t) = / {[M -M r (M)}Z(t-r{M)) + Mqz(M)} 

J M r (t) 



*(t-r(M))^dM 



(A-6) 



where qz(M) is in our case the oxygen yield for a star of 
mass M. This equation involves Z(t) inside the integral, so 
we end up with two integro-differential coupled equations. 
To solve them we run an iterative procedure assuming an 
initial value for Z(t), namely Z(t) = 0, obtaining an expres- 
sion for ez(t), and recomputing the value of Z(t) until it 
converges. Since both M r (M) and qz(M) are specified in 



the yield tables described in section 12.31 and do not have a 
simple analytical expression, all the procedure has to be run 
numerically, using interpolating functions for both functions 
and a suitable integrator. 

Once Z(t) and its inverse tz(Z) are known, we compute 
the final MDF as 



(A-7) 



MDF(log 10 (Z)) = H/(tz(Z))[l-E(ti-tz(Z))] 
dtz(Z) dZ 
dZ d(log 10 (Z)) ■ 

This is simply the aforementioned star formation 
expressed as a function of metallicity ^(tziZ)), with 
suitable variable exchange terms, and a correction term 
[1 — E(ti — tz(Z))] similar to that already introduced in 
Eqs. (|4l 1120 . This term accounts for the mass fraction of 
stars formed with a given metallicity that died before the 
final time limit t\ (i.e. the most massive ones). 
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Table 1. DRR c = {ch&, cc ■■■,CF e ) coefficients as a function of temperature 
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A composition- dependent SPH model for chemical evolution and cooling 

Table 2. Cooling function logio(Ajv) as a function of temperature T and the metallicity parameter £ 



logio(T) C=0 C=0-01 C=0-02 C=0-03 £=0.04 £=0.1 £=0.2 £=0.4 £=0.6 £=0.8 £=1 



4.0 


21.22 


21.50 


21.66 


21.78 


21.88 


22.21 


22.48 


22.75 


22.90 


22.99 


23.07 


4.1 


23.80 


23.80 


23.80 


23.81 


23.81 


23.84 


23.88 


23.95 


24.00 


24.04 


24.08 


4.2 


25.21 


25.21 


25.21 


25.21 


25.21 


25.22 


25.23 


25.25 


25.26 


25.28 


25.29 


4.3 


25.30 


25.30 


25.31 


25.31 


25.32 


25.35 


25.39 


25.45 


25.50 


25.54 


25.57 


4.4 


25.09 


25.11 


25.13 


25.14 


25.16 


25.24 


25.35 


25.50 


25.60 


25.67 


25.72 


4.5 


24.93 


24.98 


25.03 


25.07 


25.10 


25.27 


25.46 


25.67 


25.81 


25.90 


25.97 


4.6 
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